我是如何做生统作业的-2015年试题

原创 2017-06-09 徐洲更 生信媛

第一题

Lung Cancer is hard to be detected. Thus a biomarker, if successfully selected, would greatly facilitate the subsequent treatment. After primary research, a group had identified a certain candidate called protein A and they were focusing on whether the protein was a predictor of patients’ prognostic. They gathered data of patients and some healthy people from three hospitals, including the expression level of protein A together with the expression level of β-actin—-a house keeping gene, their survival time after diagnosis and others. See “instruction.docx” for details. Data is in file “lung_cancer.txt”.

Conduct a one way ANOVA to examine whether the source of data (hospital) could have impact on the expression level of protein A superficially.
Some research acclaims that whether the expression level of protein is abnormal do not have any relation with status. Based on data from hospital Alpha, could we refuse the conclusion? (Hint: Fisher’s exact test)

这道题目是不是看起来很眼熟（如果你看过我是如何做生统作业的-2014）, 这就是2014级的第四题哦，只不过第二小题不是用logistic回归，改用Fisher’s exact test而已。

注意：我这里都省略了假设部分，做题目的时候一定要写，千万要写，必须要写。

第1小题：我发现在所给真题中有去年大神提供的代码，代码部分和我之前的不同在于，注意看：

对proteinA用actin进行标准化。我觉得实际的科研中很有必要，应该充分利用数据，但是不知道答案有没有要求
在进行ANOVA分析中之前做了，正态检验，方差齐性检验。这个我没有考虑到。

注意看代码

# 设置工作环境
setwd("c:/Users/Xu/Desktop/历年生统题目/2015生统考试/")
# 导入数据
lung_data <- read.table("lung_cancer.txt", header = TRUE)
# 浏览数据
head(lung_data)
# 标准化蛋白A，这样的标准化是不是有点简单粗暴？
lung_data$A.adujst <- lung_data$exp_A /lung_data$exp_actin
# 之前在aov中因子化hospital，这里提前做
lung_data$hospital <- as.factor(lung_data$hospital)
# 看下结果
head(lung_data)
# 正态性检验
shapiro.test(lung_data$A.adujst[lung_data$hospital == 1])
shapiro.test(lung_data$A.adujst[lung_data$hospital == 2])
shapiro.test(lung_data$A.adujst[lung_data$hospital == 3])
# 方差齐性检验
bartlett.test(lung_data$A.adujst ~ lung_data$hospital)

这里节约篇幅不放正态性检验和方差齐性检验的输出，结果p值都是大于0.05的，所以可以认为数据都是正态的，不同来源的数据是等方差的。然后做方差分析，结果的p值比我大。
注意：请大家思考下到底要不要这些检验，我觉得是有必要的。但是标准化这一步，我赞同，但是助教说不做，大家按照我2014年题目来，这里提供一种思路。

exp_hospital.aov<-aov(A.adujst  ~hospital, data = lung_data)
summary(exp_hospital.aov)
#             Df Sum Sq Mean Sq F value   Pr(>F)    
# hospital      2  0.300 0.15023   7.284 0.000888 ***
# Residuals   197  4.063 0.02063                     
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

上次我到这里就结束了，上年的大神还做了多重比较, 发现第一个医院有问题。

TukeyHSD(exp_hospital.aov)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = A.adujst ~ hospital, data = lung_data)

$hospital
            diff         lwr         upr     p adj
2-1 -0.076735156 -0.13548034 -0.01798997 0.0065610
3-1 -0.078288694 -0.13703388 -0.01954351 0.0053863
3-2 -0.001553538 -0.06938663  0.06627955 0.9983886

第2小题：很简单，题目都说了fisher.test，翻到R语言实战143-144页。

fisher.test(lung_data$status[lung_data$hospital == 1], lung_data$abnormal[lung_data$hospital == 1] )

#    Fisher's Exact Test for Count Data

# data:  lung_data$status[lung_data$hospital == 1] and lung_data$abnormal[lung_data$hospital == 1]
# p-value = 0.008158
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
# 1.283453 9.026305
# sample estimates:
# odds ratio 
#   3.338326

小提示: 你知道如何选择特定行的数据么？

第二题

Researchers want to analyze the detail symptoms of Acute Keshan disease in a place. The result of phosphorus content (mg%) in patients and healthy people is showed in “Keshan_disease.txt”, in which column 1 represents phosphorus content of patients and 2 represents that of healthy people. We want to know if phosphorus content (mg%) in patients and healthy people is significantly different.

(Please record your answer, R code, and R result)

Before doing any test, please use some descriptive statistical methods to describe the difference of the data from two groups. For example, you can calculate the sample mean of the data, and draw two boxplots to show the difference of the data.
Assume that both two sample groups are drawn from normal distribution, test if the mean of the two groups are different. Note that you should test if the variances of the two groups are the same or not.
If we don’t give the assumption in b), use non-parametrical method to do the test again. Hint: Wilcoxon rank sum test or binomial test.

同志们呀，这又是一道2014年的题目呀！那年的第一题呀，这说明什么问题？这说明做历年的题目，复习平时的作业，如果感觉历年的知识点都会了，期末考随便过呀。下面的内容是我照抄我之前的答案。不对，有些地方改动了，你知道是哪里么？

明确一下关键字：描述性统计分析，方差检验，t检验，非参数检验。
请拿出《R语言实战》第二版，跟我一起翻书一起做。

第1小题请翻到第7章：基本统计分析。
首先是用summary了解他们的均值和四分位数

Keshan_disease_data <- read.csv("Keshan_disease.csv")
summary(Keshan_disease_data)
#    patient         healthy     
# Min.   :1.600   Min.   :1.280  
# 1st Qu.:4.228   1st Qu.:2.045  
# Median :4.730   Median :3.675  
# Mean   :4.652   Mean   :3.355  
# 3rd Qu.:5.385   3rd Qu.:4.228  
# Max.   :6.530   Max.   :5.780

我还需要知道一下他们的标准差

apply(Keshan_disease_data,2,sd)
# patient  healthy 
#1.058033 1.298421

各种统计函数，数学函数请翻到R语言实战的高级数据管理。这里划重点，我觉得apply是必考的内容，具体说明这里给出链接，不懂好好去看，看了还是觉得不懂的，你只能请我喝咖啡了。

画个箱图超级简单,请翻到R语言实战的基本图形：

boxplot(Keshan_disease_data)

第2小题，首先要求我们检验方差是否相同, 用到是var.test。划重点，一定要先写你的假设，不然都不叫假设检验了。
H0: 健康的人样本的方差和病人样本的方差相同
Ha: 两组的方差不同

var.test(Keshan_disease_data$patient, Keshan_disease_data$healthy)

#    F test to compare two variances

#data:  Keshan_disease_data$patient and Keshan_disease_data$healthy
#F = 0.664, num df = 39, denom df = 39, p-value = 0.2055
#alternative hypothesis: true ratio of variances is not equal to 1
#95 percent confidence interval:
# 0.3511884 1.2554349
#sample estimates:
#ratio of variances 
#         0.6639987

由于P值> 0.05,我们有理由接受原假设(听了统计的文献课，我都不知道如何看待P值了)，因此方差相同。

然后我们就可以做t.test。不知道t.test怎么用怎么办，用help(test)，

Default S3 method:
t.test(x, y = NULL,
alternative = c(“two.sided”, “less”, “greater”),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95, …)

我们需要根据自己的要求修改参数，判断是否有差异，所以 alternative = c(“two.sided”), 方差相等，所以var.equal = TRUE。
还是记住，先把假设写了
H0: 两组的均值没有差异
Ha: 两组的均值存在差异

t.test(Keshan_disease_data$patient, Keshan_disease_data$healthy, 
       alternative = c("two.sided"), var.equal = TRUE)

# Two Sample t-test

# data:  Keshan_disease_data$patient and Keshan_disease_data$healthy
# t = 4.8947, df = 78, p-value = 5.2e-06
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
 # 0.7690202 1.8234798
# sample estimates:
# mean of x mean of y 
  # 4.65150   3.35525

由于p=5.2e-06 小于0.01，所以我们有把握拒绝原假设，认为两组的均值有差异（还是不知道如何看待P值）。

第3小题的非参数检验，建议用wilcoxon rank sum或 binomial test，那我就用wilcoxon rank sum吧。请翻到基本统计分析的组间差异的非参数检验，方法就是wilcox.test，不知道怎么用，请用help(wilcox.test).
注意：写下你的假设，我这里就不写了。

wilcox.test(Keshan_disease_data$patient, Keshan_disease_data$healthy, 
            alternative = c("two.sided"))

    # Wilcoxon rank sum test with continuity correction

# data:  Keshan_disease_data$patient and Keshan_disease_data$healthy
# W = 1252.5, p-value = 1.358e-05
# alternative hypothesis: true location shift is not equal to 0

# Warning message:
# In wilcox.test.default(Keshan_disease_data$patient, Keshan_disease_data$healthy,  :
  # 无法精確計算带连结的p值

发现了一个warning message，你是不是有点慌，不要担心，只是警告而已，又不是error。如果不想看到warning添加一个exact = FALSE。

第三题：

CD38 expression is an important prognostic marker in CLL with high levels of CD38 associated with shorter overall survival. In this study, we used gene expression profiling and protein analysis of highly purified cell-sorted CD38+ and CD38- chronic lymphocytic leukemia cells to elucidate a molecular basis for the association between CD38 expression and inferior clinical outcome.

Paired CD38+ and CD38- CLL cells derived from the same patient were used to perform the analysis.
(Please record your answer, R code, and R result)
Data: GDS2676.txt GDS2676_sample.txt
Please answer the following questions:

Read in data and do normalization. Draw a boxplot using all samples before and after normalization. (Hint: limma package, function ”normalizeQuantiles”)
Find differentially expressed genes (down-regulated, CD38- < CD38+) between CD38+ and CD38- disease samples, and provide top 20 down-regulated genes. (Hint: n is small, so please use non-parameter test)
Usually you are interested in the function indicated by differentially expressed genes, for which GO enrichment is a widely used method. In order to find whether the differentially expression genes (downregulated, p<=0.05) are enriched in “leukocyte activation during immune response” (GO term), please show a conclusive result using fisher exact test. Known genes annotated with this GO are listed in “GO_2_2.txt”

这是一道新题。题目说要我们用limma, 这是一个bioconductor包，不是基础R包。这就是给我们线索了，一旦给了非基础包，我们就可以知道会出哪一类的题目了。那么如何安装limma呢。

# install limma
source("https://bioconductor.org/biocLite.R")
bioClite("limma")

我们先来学习一下limma

limmaUsersGuide(view=TRUE)

然后你会看到146页的文档，你会惊呆了的。不过不用怕，我们就看quick start，其他部分都可以略过。大体上分为以下几步：

读取数据
预处理数据，比如说标准化
拟合数据
查看结果

但是：根据题目的要求，我们连看都不用看了！！！！

基本上差异表达分析的数据读取之后都应该是如下这个样子

\	样本A	样本B	…
基因A	x	x	…
…	…	…	…

通过打开原文件，我发现体、提供的数据完全符合这个格式，根本不需要用到limma包专门的读取函数，直接用read.table就行了

第1小题：

data <- read.table("GDS2676.txt")
head(data)
boxplot(data)

我们发现有很多离群点，而且这些离群点特别不一致，这就是为什么要做标准化的原因。题目说要用normalizeQuantiles, 我们可以看下这个函数的用法

library("limma")
?normalizeQuantiles
# Usage： normalizeQuantiles(A, ties=TRUE)
# A是一个matrix， 我们导入的数据符合    
# ties：如果是TRUE，那么A的每一列都会仔细对待。好吧那就TRUE。

然后标准化，画箱图：

library("limma")
?normalizeQuantiles
data.norm <- normalizeQuantiles(data)
boxplot(data.norm)

大家的离群点特别的一致了，说明有些基因的确表达量很高。

第2小题：这时候要找差异表达的基因了，而且要用非参数的方法，而且还没有说用limma包，那就是我2014年题目中写的wilcoxon秩和检验的函数咯。

根据sample提供的数据，我们可以知道CDS-是1,3,5,7…. CDS+ 是2,4,6,8…，然后是CDS- < CDS+。还有他们来自于同一个细胞，所以应该是配对的。
这次换一种风格，用匿名函数和实体函数两种形式：

匿名函数

 p.value<-apply(data2,1,function(x)wilcox.test(x[seq(1,11,2)],x[seq(2,12,2)],paired=TRUE,alternative = "less",exact=FALSE)$p.value)

实体函数：

my.wilcox.test <- function(x){
  x <- t(x)
  CDS_sub <- as.numeric(x[seq(1,11,2)])
  CDS_add <- as.numeric(x[seq(2,12,2)])
  fit <- wilcox.test(CDS_sub, CDS_add, paired = TRUE， ternative = "less", exact = FALSE )
  pv <- fit$p.value
  return(pv)
}
pvalues.2 <- apply(data.norm, 1, my.wilcox.test)

看看结果是否一致

all(pvalues == pvalues.2)

应该返回TRUE，返回FALSE的同学，请看下参数。

然后给出TOP 20下调基因，我认为这个TOP应该是要排序的，不然直接说给出前20个结果。

sort(pvalues)[1:20]

第3小题：这一题是要做富集分析，请大家看下我写的一篇关于富集分析的文章，传送门，最后理解为什么读书会让你富有概率提高。一般而言，做GO富集分析，你要有下面这个表

\	你的结果	全部基因
注释	A	C
无注释	B	D

我们的目标是计算A,B,C,D。
首先读取GO注释：

GO_term <- read.table("GO_2_2.txt")
head(GO_term)
GO <- GO_term$V

然后提取p值小于0.05的基因。

DG.GENE <- rownames(data.norm)[pvalues.2 < 0.05]
length(DG.GENE)
[1] 617

这里说明A+B = 617

然后看DG.GENE有多少个基因在GO注释里面

 table(DG.GENE %in% GO)

FALSE  TRUE 
  607    10

所以A=10， B=607.

接着计算C，D，我们把所有用于差异分析的基因当做背景

length(pvalues)
14082

也就是说C+D = 14082，然后看着14082个基因有多少个在注释里面，有多少个不在注释里面。

table(names(pvalues.2) %in% GO)
# FALSE  TRUE 
# 13902   180

所以C=180， D=13902，那么之前的表格就是

\	你的结果	全部基因
注释	10	180
无注释	607	13902

构建fisher.test所需的矩阵。

fisher <- matrix(c(10,607, 180, 13902), nrow=2)
fisher
# [,1]  [,2]
# [1,]   10   180
# [2,]  607 13902

格式必须和上面的一模一样，不能出错。最后放上结果

 fisher.test(fisher)

#     Fisher's Exact Test for Count Data

# data:  fisher
# p-value = 0.4625
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
#  0.5967653 2.4110748
# sample estimates:
# odds ratio 
#   1.272353

从p值可以知道，并没有显著性富集呀。

limma包就用到了标准化的函数，我还以为要用里面的贝叶斯拟合呢，很尴尬，很无奈。老师果然没有高估我们。

第四题

袋中有白球5只，黑球6只，
(1) 每次无放回地抽取一个球，取三次后，这三球为两黑一白的概率；
(2) 每次有放回地抽取一个球，取三次后，这三球为两黑一白的概率；
(3) 每次无放回地抽取一个球，取三次后，这三球的顺序为黑白黑的概率；
(4) 每次有放回地抽取一个球，取三次后，这三球的顺序为黑白黑的概率；

果然是2014级的题目太难了，不然怎么会出现这个高中级别的题目。我们生统第一次课就是这个内容，所以大家记得好好回去看作业，要理解有放回和无返回的区别哦
这里我直接把之前的大神的答案搬上来了

第五题

对以往数据分析结果表明，当机器调整得良好时，产品的合格率为90%，而当机器发生某一故障时，其合格率为30%。每天早上机器开动时，机器调整良好的概率为75%，试求已知某日早上第一件产品是合格品时，机器调整良好的概率是多少？

这就是贝叶斯理论，是第2~3次课ppt的内容，这里放上百度到的定律写法

还是放上现成的答案：

感想

这次5道题目，有两道题目来自于陈老师，感觉难度下了一半。第三题目是唯一比较难的，因为用到了limma这个不太熟悉的包，不过题目说了用那些函数，压力会稍微小一点，就考验你看文档的水平如何了。然后数据管理能力还是非常重要，必须掌握如何选取数据，如何剔除数据。

最后感谢各位小伙伴的打赏支持，但是我继续不要face的放上我的二维码了。

之后会根据平时作业整理一些常用函数，不过我先休息一下吧，用脑过度，需要休息。

”FAN某”的离婚财产分割判决书（全文）

”FAN某”的离婚财产分割判决书（全文）

公益慈善｜“翼行天下一生守护”慈善项目捐赠仪式圆满举行！

哈里斯女粉搞4B运动、毒杀丈夫，回旋镖能否让美国“血流成河”

比国产光刻机更重要的IPO要来了！

我是如何做生统作业的-2015年试题

第一题

第二题

第三题：

第四题

第五题

感想

您可能也对以下帖子感兴趣

”FAN某”的离婚财产分割判决书（全文）

”FAN某”的离婚财产分割判决书（全文）

公益慈善｜“翼行天下 一生守护”慈善项目捐赠仪式圆满举行！

哈里斯女粉搞4B运动、毒杀丈夫，回旋镖能否让美国“血流成河”

比国产光刻机更重要的IPO要来了！

生成图片，分享到微信朋友圈

我是如何做生统作业的-2015年试题

第一题

第二题

第三题：

第四题

第五题

感想

您可能也对以下帖子感兴趣

公益慈善｜“翼行天下一生守护”慈善项目捐赠仪式圆满举行！